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INTERACTION OF MASSIVE BLACK HOLE BINARIES WITH THEIR STELLAR ENVIRONMENT: II. 
LOSS-CONE DEPLETION AND BINARY ORBITAL DECAY 
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ABSTRACT 



PlERO MADAU 2,3 



We study the long-term evolution of massive black hole binaries (MBHBs) at the centers of galaxies 
using detailed scattering experiments to solve the full three-body problem. Ambient stars drawn from 
a isotropic Maxwellian distribution unbound to the binary are ejected by the gravitational slingshot. 
We construct a minimal, hybrid model for the depletion of the loss cone and the orbital decay of 
the binary, and show that secondary slingshots - stars returning on small impact parameter orbits 
to have a second super-elastic scattering with the MBHB - may considerably help the shrinking of 
the pair in the case of large binary mass ratios. In the absence of loss-cone refilling by two-body 
relaxation or other processes, the mass ejected before the stalling of a MBHB is half the binary 

reduced mass. About 50% of the ejected stars are expelled in a "burst" lasting ~ 10 4 yr M 1 ^ 4 , where 
Mq is the binary mass in units of 10 6 M Q . The loss cone is completely emptied in a few bulge crossing 

timescales, ~ 10 7 yr Mg Even in the absence of two-body relaxation or gas dynamical processes, 
unequal mass and/or eccentric binaries with Mg > 0.1 can shrink to the gravitational wave emission 
regime in less than a Hubble time, and are therefore "safe" targets for the planned Laser Interferometer 
Space Antenna (LISA). 

Subject headings: black hole physics - methods: numerical - stellar dynamics 



1. INTRODUCTION 

It is now widely accepted that the formation and evo- 
lution of galaxies and massive black holes (MBHs) are 
strongly linked: MBHs are ubiquitous in the nuclei of 
nearby galaxies, and a tight correlation is observed be- 
tween hole mass and the stellar mass of the surrounding 
spheroid or bulge (e.g. Magorrian et al. 1998; Gebhardt 
et al. 2000; Ferrarese & Merritt 2000; Haring & Rix 
2004). If MBHs were also common in the past (as implied 
by the notion that distant galaxies harbor active nuclei 
for a short period of their life), and if their host galaxies 
experience multiple mergers during their lifetime, as dic- 
tated by cold dark matter (CDM) hierarchical cosmolo- 
gies, then close MBH binaries (MBHBs) will inevitably 
form in large numbers during cosmic history (Begelman, 
Blandford, & Rees 1980). Observations with the Chan- 
dra satellite have indeed revealed two active MBHs in 
the nucleus of NGC 6240 (Komossa et al. 2003), and a 
MBHB is inferred in the radio core of 3C 66B (Sudou 
et al. 2003). The VLB A discovery in the radio galaxy 
0402+379 of a MBHB system with a projected separation 
of just 7.3 pc has recently been reported by Rodriguez 
et al. (2006). MBH pairs that are able to coalesce in less 
than a Hubble time will give origin to the loudest grav- 
itational wave (GW) events in the universe. In particu- 
lar, a low-frequency space interferometer like the planned 
Laser Interferometer Space Antenna (LISA) is expected 
to have the sensitivity to detect nearly all MBHBs in the 
mass range 10 4 — 10 7 M that happen to merge at any 
redshift during the mission operation phase (Sesana et 
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al. 2005). The coalescence rate of such " LISA MBHBs" 
depends, however, on the efficiency with which stellar 
and gas dynamical processes can drive wide pairs to the 
GW emission stage. 

Following the merger of two halo+MBH systems of 
comparable mass ("major mergers"), it is understood 
that dynamical friction will drag in the satellite halo (and 
its MBH) toward the center of the more massive progen- 
itor (see, e.g., Kazantzidis et al. 2005): this will lead to 
the formation of a bound MBH binary in the violently 
relaxed core of the newly merged stellar system. As the 
binary separation decays, the effectiveness of dynami- 
cal friction slowly declines because distant stars perturb 
the binary's center of mass but not its semi-major axis 
(Begelman et al. 1980). The bound pair then hardens 
by capturing stars passing in its immediate vicinity and 
ejecting them at much higher velocities (gravitational 
slingshot). It is this phase that is considered the bot- 
tleneck of a MBHB's path to coalescence, as there is a 
finite supply of stars on intersecting orbits and the bi- 
nary may "hung up" before the back-reaction from GW 
emission becomes important. This has become known 
as the "final parsec problem" (Milosavljevic & Merritt 
2003, hereafter MM03). 

While the final approach to coalescence of binary 
MBHs is still not well understood, several computational 
tools have been developed to tackle the problem at hand. 
The orbital decay rate depends on several parameters of 
the guest binary (mass, mass ratio, orbital separation, 
and eccentricity), and on the stellar distribution func- 
tion of the host galaxy bulge. In early treatments (e.g. 
Mikkola & Valtonen 1992; Quinlan 1996, hereafter Q96), 
the stellar ejection rate and the rate of change of the 
binary semi-major axis and eccentricity were derived via 
three-body scattering experiments in a fixed stellar back- 
ground. The assumption of a fixed background breaks 
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down once the binary has ejected most of the stars on 
intersecting orbits, and the extraction of energy and an- 
gular momentum from the binary can continue only if 
new stars can diffuse into low-angular momentum or- 
bits (refilling the binary's phase-space "loss cone"), or 
via gas processes (Escala et al. 2004; Dotti, Colpi, & 
Haardt 2006). Hybrid approaches in which the rate co- 
efficients derived from numerical experiments in a fixed 
background are coupled with a model for loss-cone repop- 
ulation have been used, e.g., by Yu (2002) and MM03, 
while the limiting case in which the loss cone is constantly 
refilled but the central stellar density decreases due to 
mass ejection has been studied in a cosmological con- 
text by Volonteri, Haardt, & Madau (2003) and Volon- 
teri, Madau, & Haardt (2003). A fully self-consistent, N- 
body approach to the evolution of MBHBs, while clearly 
desirable, is limited today to N < 10 6 particles, cor- 
responding to a mass resolution of m*/M <~ 10~ 3 (e.g. 
Quinlan & Hernquist 1997; Milosavljevic & Merritt 2001; 
Hemsendorf, Sigurdsson, & Spurzcm 2002; Aarseth 2003; 
Chatterjee, Hernquist, & Loeb 2003; Makino & Funato 
2004; Bcrczik, Merritt, & Spurzem 2005). Such per- 
formance figures are not sufficient to reproduce central 
bulges, even of faint galaxies; and the small particle num- 
bers cause an artificial enhancement of star-star scatter- 
ings and of the Brownian motion of the binary, leading 
to a spurious refilling of the loss cone. 

This is the second paper in a series aimed at a detailed 
study of the interaction of MBHBs with their stellar en- 
vironment. In Sesana, Haardt & Madau (2006, here- 
after Paper I), three-body scattering experiments were 
performed to study the ejection of hypervelocity stars 
(HVSs) by MBHBs in a fixed stellar background. In 
this paper, we use a hybrid approach to investigate the 
orbital decay and shrinking of MBHBs in time-evolving 
stellar cusps. Numerically derived rates of stellar ejec- 
tions stars are coupled to an extension of the analytical 
formulation of loss-cone dynamics given by MM03. This 
method allows us to simultaneously follow the orbital de- 
cay of the pair as well as the time evolution of the stellar 
distribution function. MBHBs are embedded in the deep 
potential wells of galaxy bulges, so when the binary first 
becomes "hard" only a few stars acquire a kick velocity 
large enough to escape the host. The bulge behaves as 
a collisionless system, and many ejected stars will return 
to the central region on nearly unperturbed, small im- 
pact parameter orbits, and will undergo a second super- 
elastic scattering with the binary, as first discussed by 
MM03. Under the assumption of a spherical potential, 
we quantify the role of these "secondary slingshots" in 
determining the hardening of the pair. The plan of the 
paper is as follows. In § 2 we describe our hybrid model 
for the orbital evolution of MBHBs in a time evolving 
stellar density profile. The shrinking and coalescence of 
the binary is discussed in § 3. 

2. HARDENING IN A TIME-EVOLVING BACKGROUND 

2.1. Scattering experiments 

Our hybrid method relies on the large number of out- 
puts from the suite of three-body scattering experiments 
presented in Paper I. In the following, we briefly sum- 
marize the basic theory. Consider a binary of mass 
M = Mi + M 2 = Mi(l + q) (M 2 < Mi), reduced mass 
p = MiM 2 /M, and semimajor axis a, orbiting in a back- 



ground of stars of mass m*. In the case of a light intruder 
with to* <C M2, the problem is greatly simplified by set- 
ting the center of mass of the binary at rest at the origin 
of the coordinate system. It is then convenient to de- 
fine an approximate dimensionless energy change C and 
angular momentum change B in a single binary-star in- 
teraction as (Hills 1983) 
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Here AE/E is the fractional increase (decrease if neg- 
ative) in the orbital specific binding energy E = 
—GM/ (2a), AL Z /L Z is the fractional change in orbital 
specific angular momentum L z — \J GMa(l — e 2 ), while 
AE* and AL Z * are the corresponding changes for the in- 
teracting star. The quantities B and C are of order unity 
and can be derived by three-body scattering experiments 
that treat the star-binary encounters one at a time (Hut 
& Bahcall 1983; Q96). For each encounter one solves nine 
coupled, second-order, differential equations supplied by 
18 initial conditions. The initial conditions define a point 
in a nine-dimensional parameter space represented by the 
mass ratio q = M 2 /Mi of the binary, its eccentricity e, 
the mass of the incoming field star, its asymptotic ini- 
tial speed v, its impact parameter at infinity b, and four 
angles describing the initial direction of the impact, its 
initial orientation, and the initial binary phase. A sig- 
nificant star-binary energy exchange (i.e. characterized 
by a dimensionless energy change C > 1) occurs only for 
v < V C ^M 2 /M, where V c = ^jGMja is the binary or- 
bital velocity (the relative velocity of the two holes if the 
binary is circular, see e.g. Saslaw, Valtoncn, & Aarseth 
1974; Mikkola & Valtonen 1992). 

A set of 24 scattering experiments was performed for 
different binary mass ratios and initial eccentricities, 
each run tracking the orbital evolution of 4 x 10 6 stars. 
The binary evolution in an isotropic stellar background 
of density p and one-dimensional velocity dispersion a at 
infinity is determined by three dimensionless quantities 
(Q96): the hardening rate 
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the mass ejection rate (M e j is the stellar mass ejected by 
the binary) 

1 n ; 

(4) 
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and the eccentricity growth rate 

de 
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The hardening rate H is approximatively constant for 
separations smaller than the "hardening radius" , 

GM 2 
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(Q96). The binary is assumed to be embedded in a bulge 
of mass Mb, radius Rb, and stellar density profile ap- 
proximated by a singular isothermal sphere (SIS). Stars 
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Fig. 1. — Velocity diagram of scattered stars in three different 
speed ranges: 3V C < V < 3.2V C (black vectors), 4V C < V < 4.5V C 
(green vectors), and V > 5V C (red vectors). Upper panel: longitude 
diagram of scattered stars. Each vector length is proportional to 
the modulus of the star's total velocity (not to the velocity pro- 
jected into the xy plane). The blue ellipse shows the counter- 
clockwise orbit of the lighter black hole of the binary. Lower panel: 
latitude diagram of scattered stars. 

are counted as "ejected" from the bulge if, after three- 
body scattering, their velocity V far away from the bi- 
nary is greater than the escape velocity from the ra- 
dius of influence of the binary, rj n f = GM/ (2a 2 ). The 
SIS potential is <p(r) = -2a 2 [ln(GM B /2a 2 r) + 1] (for 
r < R B — GMb/2o ), and the escape speed from r; n f is 
then 



v csc = \/ -2<j)(r in {) 



= 2ayJ[\n(M B /M) + 1] = 5.5a, 



(7) 



where the second equality comes from the adopted bulge- 
black hole mass relation M = 0.0014 M B (Haring & Rix 
2004). Stars that do not acquire a kick velocity large 
enough to escape the host bulge, i.e. with 



V < v TCt = ^24>{R B ) - 2cj)(r in{ ) 

= 2a^f\n(M B /M) = 5.1a, 



(8) 



are allowed multiple interactions with the binary. They 
can return to the central regions on small impact pa- 
rameter orbits and undergo a second super-elastic scat- 
tering ("secondary slingshots"). Secondary scatterings 
are not allowed for stars ejected with 5.1a < V < 5.5 a, 
since even a small deviation from sphericity of the galaxy 
gravitational potential would make them miss the shrink- 
ing MBHB on their return to the center. Note that, 
even if they were able to undergo another interaction 
with the binary, such stars would not contribute sig- 
nificantly to binary hardening as long as the condition 
V > V c y/ M 2 /M — 2a -J ah /a is satisfied, since in this 
case the star-binary energy exchange would be negligi- 
ble. 

As shown in Figure [1] and discussed in details in Pa- 
per I, three-body interactions create a subpopulation of 
HVSs on nearly radial orbits, with a spatial distribution 
that is initially highly flattened in the inspiral plane of 
the binary, but becomes more isotropic with decreasing 



binary separation. The degree of anisotropy is smaller 
for unequal mass binaries and larger for stars with higher 
kick velocities. Eccentric MBHBs produce a more promi- 
nent tail of high-velocity stars and break axisymmetry, 
ejecting HVSs along a broad jet perpendicular to the 
semimajor axis. The jet two-sidedness decreases with in- 
creasing binary mass ratio, while the jet opening-angle 
increases with decreasing kick velocity and orbital sepa- 
ration. 

2.2. Loss-cone time evolution 

In the absence of loss-cone refilling by two-body re- 
laxation or other processes, the supply of stars that can 
interact with the black hole pair is limited. Analytic ex- 
pressions for non-equilibrium loss-cone dynamics based 
on the evolution of the stellar distribution function as a 
result of repeated ejections have been given in MM03. 
Here we adopt a hybrid approach instead, combining the 
results of scattering experiments with an extension of 
MM03's study. 

2.2.1. Stellar content 

The stellar content of the loss cone can be esti- 
mated from simple geometrical considerations. When 
the MBHB separation is a < ah, only a small fraction 
of bulge stars have low-angular momentum trajectories 
with pericenter distance r p < a. In a spherical velocity 
distribution, the fraction of trajectories originating at r 
and crossing a sphere of radius a < r around the center 
is 



0(r) 
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The stellar mass within the geometrical loss cone is then 

/•a i-Rb 

M* = / 4?rr 2 p(r) dr + / 4?rr 2 p(r)6(r) dr. (10) 

JO J a 

For a SIS p(r) — a 2 /(2irGr 2 ), and equation (fT0|) is read- 
ily integrated to yield in the limit a <C R B 

7T(T 2 7T / a \ „ 

M* ~ — o= - — M 2 . 
G 4 \a h J 

The above scheme is oversimplified, as it assumes stel- 
lar trajectories to be straight lines. The gravitational 
field of the stellar mass distribution increases the net 
number of distant stars with pericenter distances r p < a. 
Consider a star at distance r > a from the binary mov- 
ing with random velocity v. For an SIS, conservation of 
energy gives: 



(11) 



v 



,„ ,-4a 2 ln(r», (12) 

where v p is the star's velocity at pericenter. If b is the 
impact parameter at distance r, angular momentum con- 
servation yields 



r 2 p [l + 4(a 2 /v 2 )Hr/r p )}. 



(13) 



The second integral on the right-hand-side of equation 
(flTJl) can then be rewritten as 



47rr 2 p(r) 0(r) 



4tt v 2 f{v) [1 + 4(a 2 /i; 2 ) ln{r/a)}dv \ dr, 



(14) 
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where f(v) is the stellar velocity distribution. For a 
Maxwellian, the above equation can be simplified by set- 
ting v = (v) = \fia in equation (|13p : one can then define 
a 0-factor that includes gravitational focusing 
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(15) 

Numerical integration of equation (JTUJ) finally yields for 
the stellar mass in the loss cone 



8.2a 2 „ fa , , . 
M* ~ — — a ~ 2 — ) M 2 . 



(16) 



Note that a fraction ~ O.5M2 of the mass contained in 
the loss cone when the binary becomes hard (a = a^, 
M* ~ 2M2) lies within a^. Let t = be the time at 
which the binary separation is a = ah- The number flux 
of stars into the geometrical loss cone, i.e. the flux of 
stars with r > au at t = that interact with the binary 
at a later time t, is 



T 



2 V / 3V 3 



6, 



(17) 



where 6 = Q(^/3at + o^). 



2.2.2. Energy exchange 

Let us denote with £ the total binding energy of the 
MBHB, £ = GM 1 M 2 /2a. The total energy transfer rate 
from the binary to stars in the loss cone can be written 

as 



d£_ ^ AE*(a)M*(a) 
dt i c h 



(18) 



where t c h is a characteristic interaction timescale and 
AP*(a) ~ G/i/a is the characteristic specific energy gain 
of stars as a consequence of the gravitational slingshot. 
MM03 have written equation ([18)1 in the case of re- 
turning stars, i.e. kicked stars that do not escape the 
host bulge and can have a secondary super-elastic in- 
teraction with the MBHB. Returning stars have energy 
E(a) ~ <f>(ri n f) + AE*(a), and their interaction timescale 
t c h can be identified with the typical radial period of stars 
in an SIS potential, 

t ch ~ P(E) = P(0) exp( J B/2 ( 7 2 ) 
= P(0)(Af/Af s )exp 



2C 



(t)- 1 



(19) 



Here, P(0) = ^/ttGAIb /2a 3 is of order the bulge crossing 
time, and the average dimensionless energy change C 
is of order unity and nearly independent of a for a < 
a h . As noted by MM03, in this case M*(a)A_E*(a) oc 
a a ~const. From equation (|16j) simple calculations 
lead to 
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where ai is the binary separation at time t = t\ when 
secondary slingshots start, and Pi is the period of stars 
with energy E{a\). 




10- 3 0.01 

t/P(0) 

Fig. 2. — Decay of binary separation a (in units of a^) as 
a function of time [in units of the bulge crossing time P(0) = 

1.32 X 10 7 yr Afg /4 for a SIS]. The curves, from bottom to top, 
are for g = 1/243, 1/81, 1/27, 1/9, 1/3, 1. The q = 1 case is com- 
pared to the analytical estimate (long-dashed line) of MM03, and 
to an N-body simulation (short-dashed line) of MM03 performed 
with 18,000 stars initially in the loss cone, and the stellar potential 
replaced by a smooth component to prevent rel axa tion. We set 
M = 0.0014 Mg and use the M — a relation (ea.[24t. 

MM03's analysis can be expanded to account for the 
effect of three-body scatterings when the MBHB first 
becomes hard at separation a — a^. Stars with r < a^ 
at t = will interact with the binary within a timescale 
£ c h ~ a/ l /( v / 3cr). Substitution of t c h in equation (fT5)l. 
followed by simple algebra, leads to the expression 



(21) 



A further contribution to the shrinking is associated with 
stars having r > at t — that are bound to enter the 
loss cone at later times. This population has total mass 
~ 1.5 M2, and its contribution to the orbital decay is 
given by 

AP»(a)m*^-"7ra 2 




d£ 



(22) 



A straightforward substitution gives 

V3cr 



dt\ a ) VI 




e(V3at + a h ). (23) 



The above equation holds for a bulge crossing time, t < 
Rb/VSct, and must be solved numerically. 

2.2.3. Orbital decay 

The simple analytical formulation described above can 
be refined using results from our scattering experiments 
(Paper I) . This allows us to follow at the same time both 
the orbital decay of the binary and the evolution of the 
distribution function of interacting stars. The procedure 
is the following. We first isolate, from the initial distri- 
bution of kicked stars, the new loss cone, i.e. the subset 
of stars returning to the center on orbits with r p < a\, 
where a\ is the binary separation at the end of the first 
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TABLE 1 

Binary hardening: the impact of returning stars 



q 


ah/ai 


Oft /a j 


1 


2.81 


4.41 


1/3 


2.19 


3.07 


1/9 


1.64 


2.09 


1/27 


1.29 


1.49 


1/81 


1.12 


1.19 


1/243 


1.04 


1.06 



Note. — Binary shrinking factors, /ai is the binary shrinking 
after the first interaction only, while a^/af take into account for 
subsequent reejections up to the fourth interaction. 

interaction with the stellar background. Then we com- 
pute the hardening rate H by averaging H\(v) (provided 
by our scattering experiments) over the velocity distribu- 
tion function of such stars, which are allowed to interact 
with the binary for the timcscalc in equation (|19[) , again 
averaged over the stellar velocity distribution. After each 
step, the velocity distribution function of returning stars 
is updated, and so is the timescale of the following in- 
teraction. We iterate the process until the loss cone is 
emptied. Convergence to the final stalling separation is 
usually obtained after > 4 iterations. The mathemati- 
cal details of the numerical procedure are given in the 
Appendix. 

In order to specify the two parameters defining the 
SIS, the Haring & Rix (2004) bulge-black hole mass re- 
lation was complemented by the M — a relation (Fer- 
rarese & Merrit 2000; Gebhardt et al. 2000) proposed 
by Tremaine et al. 2002: 

(7 7 o = 0.84 M 6 1/4 , (24) 

where 070 is the stellar velocity dispersion in units of 70 
km/s, and M 6 is the MBHB mass in units of 10 6 M Q . 
Using equation (|24|) . the hardening radius can then be 
written as 

a h ~ 0.32 pc M^ 2 ^^, (25) 

where M^fi = M2/1O 6 M0. The binary separation as a 
function of time is shown in Figure [2J where it is also 
compared to the results of an N-body simulation and 
to an analytical prescription, both presented in MM03 
(their Fig. 6). The agreement with the simulation is 
fairly good. Our results, in terms of the final separation 
achieved as a function of g, are perfectly consistent with 
the stalling radii estimated by Merritt (2006). 

Figure [5] shows how the rate of orbital decay de- 
clines after a few bulge crossing times ~ P(0) = 1.32 x 

10 7 yr Mg^ 4 : this is due to the decreasing supply of low 
angular momentum stars from the outer regions of the 
bulge, once the stars in the central cusp have interacted 
with the binary. Note that equal mass binaries shrink 
more than unequal binaries: this is both because of the 
scaling of the stellar mass available for the interaction 
and of the energy exchanged during a typical three-body 
encounter. It is easy to see that d(l/a)/d\nt oc q. Ec- 
centricity plays a marginal role in the orbital evolution of 
the MBHB. For a given q, the orbital shrinking is larger 
by at most 10% for highly eccentric binaries. 

Figure [3] shows an example of the role of secondary 
slingshots on orbital shrinking for an equal mass circular 



binary. The lower panel clearly illustrates how successive 
interactions of stars returning on quasi radial orbits can 
reduce the final binary separation by an extra factor of 
order 2, i.e. a/ ~ ai/2. The progressive emptying of the 
loss cone is sketched in the upper panel, where we plot 
the (differential) mass in stars approaching the binary 
with a given periastron. After the first interaction only 
few stars are kicked out from the bulge. The loss cone, 
while substantially hotter, remains nearly full and only 
gets progressively depleted as secondary slingshots take 
place. Table [T] quantifies the role of returning stars for 
different values of the binary mass ratio q: returning 
stars can increase the shrinking of the MBHB by as much 
as a factor of 2, and play a larger role for equal mass 
binaries. This is because binary-to-star energy exchange 
is significant only for V < ^fqV c (see Paper I). After 
the first binary-star interaction, the stellar population is 
heated up and stars have on average V ~ y/qV c . In the 
case q = 1, the binary shrinks by a significant factor, V c 
increases, and most of the returning stars have V ;$ V c - 
the hardening process is still efficient. By contrast, when 
q <C 1, V c does not increase appreciably after the first 
interaction, returning stars have V ~ \foVci an d binary 
hardening stops. 

3. DISCUSSION 

Under the assumed criterion for stellar ejection, we can 
compute M e j, the mass of stars expelled with V > v csc . 
An example of the effects of the slingshot mechanism on 
the stellar population is shown for an equal mass circu- 
lar binary in Figure IH where the initial (t = 0) velocity 
distribution of interacting stars is compared to the distri- 
bution after loss-cone depletion. As already mentioned 
in the previous section, after the first interaction with the 
binary a large subset of kicked stars still lies in the loss 
cone of the shrinking binary, has velocities v < w ret , and 
is potentially available for further interactions. These 
are the stars we termed "returning" . While scattering 
with the binary increases the stellar velocity thus reduc- 
ing the energy exchanged in secondary interactions, it 
moves kicked stars on more radial orbits thus reducing 
their impact parameter as well. Our calculations show 
that the high velocity tail of the distribution depends on 
the MBHB eccentricity e, although the effect of changing 
e is small for small values of q. In this case, fewer stars 
are kicked out compared to the case q < 1, but at higher 
velocities on the average. In general, both a small mass 
ratio and a high eccentricity increase the tail of HVSs. 
Integrating the curves in Figure 3] over velocity gives the 
mass of interacting stars: this is ~ 2M for the case shown 
(q = 1), ~ 1.2M for q = 1/3, and ~ 0.6M for q = 1/27 
(all assuming e = 0; we checked that eccentricity plays a 
negligible role). 

Figure [5] depicts the ejected mass M C j normalized to M 
(left scale) and to M2 (right scale), as a function of q. Our 
results show that M ej /M ~ 0.5^/Af = 0.5g/(l + q) 2 , i.e., 
M e j/M 2 ~ 0.5/(1 + q), both ratios being independent of 
the total binary mass. The rate of stellar mass ejection is 
shown in Figure[6]as a function of time. A fraction < 50% 
of the expelled stars is ejected in a initial burst lasting 
~ a/i/c, this fraction being smaller for smaller binary 
mass ratios. The burst is associated to the ejection of 
those stars already present within the geometrical loss 
cone when the binary first becomes hard. Note that, for 
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to- 6 10- 5 10- 4 10- 3 0.01 0.1 
t/p(o) 

Fig. 3. — Upper panel: evolution of the loss-cone population in 
terms of the (differential) stellar mass that approaches the binary 
with a given periastron r p . Thin line: initial loss cone. Solid lines, 
from top to bottom: loss-cone population after the 1st, 2nd, 3rd, 
and 4th interaction. Lower panel: binary separation as a function 
of time. From bottom to top, the curves depict the shrinking as- 
sociated with only the first one, two, three, and four interactions, 
respectively. An equal mass, circular binary is assumed. 




0.01 



Fig. 4. — Stellar velocity distribution for an equal mass, circular 
binary, at different stages of binary hardening. The vertical lines 
mark v TC t = 5.1a (we recall that v eBC = 5.5cr). Dashed lines: from 
top to bottom, distribution of stars in the shrinking loss cone before 
the 1st, 2nd, 3rd and 4th iteration. For clarity, the initial loss cone 
distribution is marked with a thicker line. Thin solid lines: from 
top to bottom, distribution of stars that have received 1, 2, 3 and 
4 kicks. Thick solid line: final stellar velocity distribution after 
loss-cone depletion is completed. 

small q, mass ejection is already significant at a ~ ah, 
as in this case the binary orbital velocity is V c > v esc for 



< 



ah- 



Is the amount of ejected mass sufficient to shrink the 
MBHB orbit down to the GW-dominated regime? To 
answer this question, we start defining the "final separa- 
tion" af as the separation reached by the binary before 
complete loss-cone depletion, i.e. after a few bulge cross- 
ing time (~ 10 7 yrs, weakly depending on binary mass 



as oc M 1 ! 4 ). We must then compare a/ to the separation 
at which the orbital decay timescale from GW emission, 



5c 5 



256G 3 MiM 2 MF(e) 
MM 1 M 2 



0.25 Gyr 



10 18 - 3 M m 3 



F{e)^ 



O.OOlpc 



(26) 

(Peters 1964), is shorter than, say, 1 Gyr. Here, to 4th 
order in e, 



F(e) 



(1 



1 73 z 

1 + 24 6 + 



37 
96 ( 



(27) 



Inverting equation (|26p . one can define the separation 
oqw at which the binary will coalesce in a given time i, 



256G 3 

—tM 1 M 2 MF{e) 



1/4 



/ MMi m 2 y/' yi 

Vio 18 - 3 m V K ' 9 



(28) 



5c 
0.0014 pc 



where tg = i/lGyr. Using equations (f2"S"| and 
one finds that a//a G w oc M~ 1 / 4 q 3 / 4 (see also MM03, 
eq. 90), i.e. the more massive the binary and the smaller 
the binary mass ratio, the smaller the factor the binary 
must shrink to reach the GW emission regime. Eccen- 
tricity plays a double role. For a given binary mass and 
mass ratio, on one hand the hardening rate slightly in- 
creases with increasing eccentricity (< 20% from e = 
to e — 0.9), leading to a smaller a/; on the other hand, 
from equations ([27]) and ([28]) . «gw is larger for larger 
e thus reducing the at-acw gap- An important effect, 
included in our calculations, is that the eccentricity typi- 
cally increases during the binary-star interaction, though 
the functional form of i^e) is such that the effect is sig- 
nificant only for binaries with e > 0.6 already at a^. In 
other words, for MBHBs with initially low eccentricities, 
the increase of e during the gravitational slingshot affects 
only weakly the final a / /acw ratio. 

It is interesting to compare the total mass actually 
ejected prior to complete loss-cone depletion with the 
stellar mass that must be expelled in order to reach a fi- 
nal orbital separation where tew (a/) = 1 Gyr, i.e. where 
GW emission leads to coalescence within 1 Gyr. An ex- 
ample is shown in Figure [5] The shaded areas define 
such mass (in units of M) for a 10 6 M Q and a 10 9 M Q 
MBHB, where the top boundary assumes e = and the 
bottom e = 0.9. Note how the e = 0.9, M = 10 6 M Q 
lower boundary practically coincides with the e = 0, M = 
10 9 Mq upper one. The figure clearly shows how, even 
in the absence of other mechanism driving orbital decay, 
pairs involving genuinely supermassive holes should not 
stall, while for lighter binaries both a small mass ratio 
and a large eccentricity are probably required for coales- 
cence to take place. 

Using our hybrid model, we can also sample the 
{Mi, g, e) 3-D space, compute the separation a/ and the 
eccentricity e at a/, then fold the calculated values of 
a / and e into equation (|26p , and finally compare £gw to 
the Hubble time at two reference redshifts, z = 1 and 
z = 5. In Figure [71 binaries that will coalesce within a 
then Hubble time after loss-cone depletion populate the 
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Fig. 5. — Ejected stellar mass M c j normalized to the total binary 
mass M (left scale, solid points), and to the mass of the lighter 
binary member M2 (right scale, empty points), as a function of 
binary mass ratio. The curves are polynomial interpolations. Note 
that the ratios M e j/M and M G j/M2 do not depend on the absolute 
value of M, and are nearly independent on e. Upper, dark shaded 
area: the mass (normalized to M) a M = 10 6 Mq binary needs to 
eject to reach a final separation aj such that tow = 1 Gyr. Top 
and bottom boundaries assume e = and e = 0.9, respectively. 
Lower, light shaded area: same but for a M = 10 9 Mq binary. 
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t/P(0) 

Fig. 6. — Ejected stellar mass per unit logarithmic time interval, 
as a function of time. A circular binary is assumed. Solid line: 
q = 1. Long-dashed line: q = 1/3. Short-dashed line: q = 1/9. 
Dot-dashed line: q = 1/27. 

diagonally-shaded area in the Mi — q plane, while the 
vertically-shaded area marks MBHBs that, if driven to 
coalescence by z = 1 or z = 5, would be resolved by 
LISA with a signal-to-noise ratio S/N > 5 (see Sesana 
et al. 2005 and references therein for details). The region 
of overlap selects unequal mass, highly eccentric MBHBs 
with M > 10 5 that can shrink down to the GW emission 
regime in less than an Hubble time, and that are "safe" 
targets for LISA even in the pessimistic case, treated 
here, of stellar slingshots+loss-cone depletion with no re- 
filling. 
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Fig. 7. — Mi — q plane. The vertical shaded area shows LISA 
potential targets with S/N > 5. The diagonal shaded area on the 
lower right corner marks binaries that will coalesce within a then 
Hubble time after loss-cone depletion. In each panel, the assumed 
redshift and eccentricity of the MBHBs are labeled. 

It is important to remark, at this stage, that our cal- 
culations are meant to define a minimal model for the 
evolution of MBHBs, and that several other mechanisms 
may help the orbital decay and widen the range of poten- 
tial LISA targets. First, we have assumed all the stars in 
the loss cone to be unbound to the MBHs. In a realistic 
case, each MBH will bind stars inside its radius of influ- 
ence rinf: the star binding energy can be extracted by 
slingshot, hence enhancing binary hardening. This effect 
is not expected to be important for equal mass binaries 
as, in this case, ~ n n f , and only a small fraction of in- 
teracting stars will be bound to the binary. Indeed, our 
results match well the numerical simulations of MM03 
(Fig. [2]) . For lower mass ratios, however, it is <C nnf 
and most stars in the loss-cone are actually bound to the 
binary. A forthcoming paper will be devoted to an anal- 
ysis of three-body scattering experiments for a MBHB 
with bound stars, providing a more realistic model for 
the case q <C 1. 

Second, even in spherical stellar bulges, loss-cone refill- 
ing due to two-body relaxation (Yu 2002; MM03) and the 
wandering of the black hole pair in the nucleus (Quinlan 
& Hernquist 1997; Chatterjee et al. 2003) could both 
increase the amount of stellar mass interacting with the 
MBHB. The two-body relaxation timescale is such that 
loss-cone refilling is probably important for M < 10 6 M Q . 
The Brownian motion timescale is of the order of 15 Gyr 
for a 10 6 Mq binary, and scales as M 5 / 4 . It is then likely 
that the two effects considered here may affect orbital 
decay only for light binaries, helping them to cover the 
residual gap between a / and acw and leading light bina- 
ries to coalesce within a then Hubble time even at high 
redshifts. On the other hand, their contribution to the 
shrinking of supermassive binaries with M > 10 6 M Q is 
probably negligible. 

If the stellar bulge is not spherical, but axisymmetric, 
stars on highly eccentric orbits are typically centrophilic 
(Touma & Tremaine 1997; Magorrian & Tremaine 1999). 
In this case, the loss cone is substituted by a "loss wedge" 
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(see Yu 2002 for a detailed discussion). The stellar con- 
tent of such wedge is larger than that of the correspond- 
ing loss cone, and depends on the degree of flattening e 
of the stellar distribution. Typically, for a galaxy with 
e = 0.3, the stellar content of the loss wedge is a order of 
magnitude larger than the stellar content of the loss cone 
were the bulge spherical (Yu 2002). Note that Faber et 
al. (1997) estimate an average e = 0.36 for a sample of 
galaxies, leading to the conclusion that the hardening of 
a MBHB in such a potential might be much faster than 
our "spherical" estimate. We also recall that triaxial po- 
tentials drive many bulge stars on chaotic orbits, many 
of them centrophilic: one then expects an increase in the 
number of interacting stars similar to that produced by 
axisymmetric potentials (Merritt & Poon 2004; Berczik 
et al. 2006). 

Finally, MBHB orbital evolution can also, at least par- 
tially, be driven by drag in a gaseous nuclear disk. The 
role of gas is, basically, twofold. On 100 pc scales, the 



disk drastically increases dynamical friction, reducing the 
timescale on which MBHs can reach the center of the 
bulge (Escala et al. 2004; Dotti et al. 2006). On par- 
sec scales, torques induced by the disk can drive the bi- 
nary to decay on a timescale of order the gas accretion 
time (Ivanov et al. 1999; Armitage & Natarajan 2002). 
It is important to point out that the interaction with a 
gaseous disk typically circularizes the binary orbit (Dotti 
et al.2006), hence maximizing the a^ — oqw S a P- If this is 
the case, the slingshot driven coalescence would be more 
difficult to achieve. 
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APPENDIX 

NUMERICAL INTEGRATION OF THE MBHB ORBITAL DECAY 

The average hardening rate for a Maxwellian stellar velocity distribution f(v,a) — (27r<7 2 )~ 3 / 2 exp(— v 2 /2a 2 ) is 
given by 



where 



H{a) = / f(v,cx)-Hi{v)AW 
Jo v 



H t (v) = 8tt / (C)xdx 
Jo 



dv, 



(Al) 
(A2) 



is the dimensionless hardening rate if all stars have the same velocity v, x = b/ \/2GMa/v 2 is the dimensionless impact 
parameter, and the ene rgy exchange (C) is averaged over the orbital angular variables (Paper I; Q96). An expression 
analogous to equation (lAlj) relates the thermally-averaged eccentricity growth rate K(a) to K\{v): 



= (l-e 2 )J-(B-C)xdx 
2e f °°(C)xdx ' 



(A3) 



where (B — C) is the mean angular momentum minus energy exchange. For a binary with given mass ratio and 
eccentricity, the quantities C and B, and thus Hi and K\, are only function of the ratio v/V c cx vy/a, where V c is the 
binary circular velocity. Given an incoming velocity v, we record the bivariate distribution hi(V,b'\v) of stars with 
ejection speeds in the interval V, V + dV, leaving the binary with an "exit" impact parameter in the interval b' , b' + db' . 
The distribution function is normalized so that 



hi(V,b'\v)dVdb' 



1. 



(A4) 



The subscript "1" indicates that the scattering experiments are performed for a binary at separation a = 1. 

The interaction of the MBHB with stars in the loss cone is assumed to take place in discrete steps. The binary 
first interacts with a given population of stars and shrinks accordingly. We then isolate the returning sub-population, 
which becomes the input for the next step, and so on. Consider the binary at separation at, interacting with a stellar 
population of mass M*^ and (normalized) velocity distribution fi(v). The orbit decays according to the differential 
equation 

GP H(a), (A5) 



dt 



where 



H(a) 



< v > 



fi(v) <v> Hi(vy/a) dv, 
o v 



and p is the stellar density. Straightforward integration of equation (|A5|) gives 

< v > f a ' da' 



t(a) = 



Gp J a a' 2 H(a'Y 



(A6) 



(A7) 
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the time the orbit needs to shrink from to a. The solution above is not physically meaningful, as the time variable 
involved depends on the particular value assumed for p. Stated more directly, to solve equation (|A5j) we need to set a 
realistic pace at which star-binary interactions occurs. This can be done in two steps. 
First, we write the stellar mass M* that will interact with the binary in a given time t as 

M*(t) = k*t, (A8) 



where 



and the maximum allowed impact parameter is 



Mv)nb max (v) pvdv, (A9) 



VLJV) =«?(!+ 2 ~^f) ■ (A10) 



'J 

Equation (|A8j) is valid as long as M*(t) < M*^. Now we simply write t = Af*/fc„, and substitute into equation (|A7|) . 
which now gives M*(a), the stellar mass interacting with the binary as the orbit shrinks from aj to a. Note that in 
M*(a) term p cancels out, i.e., the interacting mass is independent of any pre-assigned value for the stellar density. 
We can now numerically invert the equation for M*(a), obtaining a(M*), i.e., the binary separation as a function 
of the interacting mass. We define the final separation a/ = a(M* ; j). The very same procedure can be applied for 
the evolution of the binary eccentricity e, resulting in a function describing e as a function of the interacting mass, 
e = e(M*). 

We need now to relate M* to physical time. The stellar mass interacting with the binary per unit time is 

dM* dM * dv „ . . , r dv , . „ „ . 

" ^ =/i(u)M *^< < All > 
where the term dv/dt can be computed considering the typical interaction time for stars with velocity v in a SIS 
poten tial, i.e., t(v) = P(0)exp[(t> 2 — u 2 ct )/4cr 2 ]. Straightforward algebra yields dv/dt as a function of t. Equation 
(|A11|) can be then integrated, and the resulting M*(t) finally substituted into a(M„) to give the time evolution of the 
binary separation a(t). 

We can now compute the distribution of stars that, after the interaction, have velocities V < v rct , and are then 
available for a further encounter with the MBHB. The bivariate distribution h a (V,b'\v) can be extracted from the 
hi(V,b'\v) distributions recorded in the scattering experiments. In the three-body problem integration, the binary 
mass and separation are taken as unity. This set-up allows to rescale the velocities and trajectories of kicked stars for 
any given physical value of the binary mass and separation. It can be shown that 

h a (V,b'\v) = x hx \^=,b'a\vy/a] , (A12) 



where the prefactor 1 / normalizes the distribution according to equation (|A4I) . The normalized distribution of 
scattered stars must be averaged over the input velocity v and the separation a, and can be written as 

C io°° ft v ^ ( dM * / dQ ) h - ^ b » dv da 
Jo°° Jo°° la! fo°° fi(a,v)(dM*/da)h a (V,b>\v)dvdadVdV ( J 

Here, the distribution of stars with velocity between v and v + dv that are going to approach the binary within a 
distance < a is 

fM= , (A14) 



b\v,a)=a'\l+'^f), (A15) 



where 

/ 

h 2 („ n\ - n 2 I 1 4- : 

av 2 J 

and fe max (w) is given by equation (|A10|) . The change of the interacting mass with binary separation dM^/da is obtained 
by differentiation of the function M^a), obtained above. 
Finally, the (normalized) velocity distribution of returning stars f r (v) can be computed as 

Io f h ^ h ') dh ' 

Sr{v) = h(v,b')db'dv ' ( A16) 

and the mass avaiable for the subsequent interaction M*^ is given by 

M* r = M* i Jo r< d rOO , ' • A17 

■ / °° J °°h(V,b')db'dV K ' 

The numerical procedure can be then iterated, considering the binary at a starting separation a t — a,f, interacting with 
a stellar population whose (normalized) velocity distribution is fi(v) = f r (v), allowing a total mass of interacting stars 
M* t i — M*^ r . For the first interaction only, we assume that the stars already within the binary separation interact on 
a time scale ~ a^/cr, along the lines discussed in Section 2.2.2. 
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